/// by Alvaro Cuno, Fev2008

#ifndef __RBF_EVAL_KERNEL_3TO3_H_
#define __RBF_EVAL_KERNEL_3TO3_H_

#include "rbfEvalKernel_3to1.h"

class RbfEvalKernel3to3 : public RbfEvalKernel3to1 {
   
public:
   
   RbfEvalKernel3to3() {
      
      RbfFunction();
      
      npolyterms = 4;
      pointsTexSize = 0;
      isDataInTexture0 = true;
      evalKernel.fragment_source("rbfEval_3to3.fsc");
      fbo.setNumberOfChannels(3);
   }
   
   virtual ~RbfEvalKernel3to3() {
      glDeleteTextures(1, &texPoints);
   }
   
   void setCentersAndWeights(vector<Point3> &c, vector<Point3> &w) {
      centers = c;
      weights = w;
   }
   
   void execute(vector<Point3> &values) {
      
      fbo.on();
      isDataInTexture0 = true;
      evalKernel.use(true);
      
      /// Transfers the points to the GPU memory
      glActiveTexture(GL_TEXTURE2);
      glBindTexture(GL_TEXTURE_RECTANGLE_NV, texPoints);
      evalKernel.set_uniform("texturePoints", 2);
      
      glActiveTexture(GL_TEXTURE0);
      for (uint i=0; i<centers.size(); ++i) {
         glEnable(GL_TEXTURE_RECTANGLE_NV);
         
         if (isDataInTexture0) {
            glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT); // writting into texture 0
            glBindTexture(GL_TEXTURE_RECTANGLE_NV, fbo.getTextureId(1)); // reading from texture 1
         }
         else {
            glDrawBuffer(GL_COLOR_ATTACHMENT1_EXT);
            glBindTexture(GL_TEXTURE_RECTANGLE_NV, fbo.getTextureId(0));
         }
         evalKernel.set_uniform("center", centers[i][0], centers[i][1], centers[i][2]);
         evalKernel.set_uniform("weight", weights[i][0], weights[i][1], weights[i][2]);
         evalKernel.set_uniform("texturePhi", 0);
         isDataInTexture0 = !isDataInTexture0;
         
         glBegin(GL_QUADS);
           glVertex2f(0, 0);
           glVertex2f(pointsTexSize, 0);
           glVertex2f(pointsTexSize, pointsTexSize);
           glVertex2f(0, pointsTexSize);
         glEnd();
         
         glDisable(GL_TEXTURE_RECTANGLE_NV); 
      }
      glFinish();
      
      isDataInTexture0 = !isDataInTexture0;
      evalKernel.use(false);
      
      downloadData((Number*)&values[0]);
      fbo.off();
   }
   
   ///
   void downloadData(float *data) {
      
      if (isDataInTexture0)      { glReadBuffer(GL_COLOR_ATTACHMENT0_EXT); }
      else                       { glReadBuffer(GL_COLOR_ATTACHMENT1_EXT); }
      
      glReadPixels(0, 0, pointsTexSize, pointsTexSize, GL_RGB, GL_FLOAT, data);
      glFinish();
   }
   
   void addPolynomialTerm(vector<Point3> &pEval) {
      
      const int &n = centers.size();
      for (uint ii=0; ii<pPoints->size(); ++ii) {
         const Point3 &p = (*pPoints)[ii];
         Number rx=0, ry=0, rz=0;
         for (unsigned int i=0; i<npolyterms; ++i) {
            const Number &valpoly = polynomial(p, i);
            rx += weights[n+i][0]*valpoly;
            ry += weights[n+i][1]*valpoly;
            rz += weights[n+i][2]*valpoly;
         }
         pEval[ii][0] += rx;
         pEval[ii][1] += ry;
         pEval[ii][2] += rz;
      }
   }
   
   /// Performs the product: w = iM*fcenters
   /// @param iM the RBF inverse matrix
   /// @param cntrs the centers array
   /// @param w the weights array
   void genWeights(vector<Point3> &fcenters) {
      
      assert(centers.size() == fcenters.size());
      unsigned int n = fcenters.size();
      assert(npolyterms == iM.m-n);
      Matrix RightSide(iM.m, 3);
      for (unsigned int i=0; i<n; ++i) {
         const Point3 &Q = fcenters[i];
         RightSide(i,0)=Q[0];
         RightSide(i,1)=Q[1];
         RightSide(i,2)=Q[2];
      }
      
      /// additional restrictions initialization
      for (unsigned int i=0; i<npolyterms; ++i) {
         RightSide(n+i,0)=0.0;
         RightSide(n+i,1)=0.0;
         RightSide(n+i,2)=0.0;
      }
      Matrix lambdas = iM*RightSide;
      
      weights.resize(iM.m);
      for (unsigned int i=0; i<weights.size(); ++i) {
         weights[i][0]=lambdas(i,0);
         weights[i][1]=lambdas(i,1);
         weights[i][2]=lambdas(i,2);
      }
   }
   
//~ protected:
   vector<Point3>    weights;       // the vector of weights
};

#endif
